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In this paper, we develop a lattice Boltzmann model for relativistic magnetohydrodynamics 
(MHD). Even though the model is derived for resistive MHD, it is shown that it is numerically 
robust even in the high conductivity (ideal MHD) limit. In order to validate the numerical method, 
test simulations are carried out for both ideal and resistive limits, namely the propagation of Alfven 
waves in the ideal MHD and the evolution of current sheets in the resistive regime, where very good 
agreement is observed comparing to the analytical results. Additionally, two-dimensional magnetic 
reconnection driven by Kelvin-Helmholtz instability is studied and the effects of different parameters 
on the reconnection rate are investigated. It is shown that the density ratio has negligible effect on 
the magnetic reconnection rate, while an increase in shear velocity decreases the reconnection rate. 
Additionally, it is found that the reconnection rate is proportional to , a being the conductivity, 
which is in agreement with the scaling law of the Sweet-Parker model. Finally, the numerical model 
is used to study the magnetic reconnection in a stellar flare. Three-dimensional simulation suggests 
that the reconnection between the background and flux rope magnetic lines in a stellar flare can 
take place as a result of a shear velocity in the photosphere. 

PACS numbers: 47.ll.-j, 47.65.-d, 47.75.-|-f 


I. INTRODUCTION 

Magnetic fields are an essential component of many as- 
trophysical phenomena, such as relativistic jets [l| , active 
galactic nuclei Q, gamma ray bursts Q, pulsar winds 0, 
and stellar flares [^. Since in most of these phenomena 
the plasma is electrically neutral and the characteristic 
times between collisions are much smaller than the typi¬ 
cal time scale of the system, the magnetohydrodynamics 
(MHD) approximation is appropriate. Due to the fact 
that relativistic effects play a major role in the dynamics 
of these phenomena, relativistic MHD description is of 
special interest in this perspective. Except for some sim¬ 
ple geometries, most of the studies are based on numer¬ 
ical simulations, since the equations of relativistic MHD 
are extremely difhcult to solve analytically. 

Ideal MHD is defined as the limit where the electri¬ 
cal conductivity a goes to infinity (electrical resistivity 
ry = l/cr —>■ 0). In this framework many numerical mod¬ 
els have been developed over the last decade dealing with 
the ideal relativistic MHD [^-Q. As it will be explained 
later, the ideal MHD assumption not only makes the solu¬ 
tion of the relativistic MHD considerably simpler, but it 
is also a fairly good approximation for many high-energy 
phenomena. However, in several situations such as neu¬ 
tron star mergers @ or central engines of gamma ray 
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bursts 0: the conductivity can be small and the ideal 
MHD assumption is not valid any longer. 

More importantly, magnetic reconnection only takes 
place when resistivity exists in the plasma. This process 
is the driver of explosive events in astrophysical plas¬ 
mas, in which magnetic field lines break and reconnect 
and the magnetic field topology goes through a sudden 
change. During this process, plasma releases the mag¬ 
netic energy and converts it into thermal and kinetic en¬ 
ergy on a short timescale. Magnetic reconnection has 
been proposed to have an influential role in many as¬ 
trophysical observations, namely as a cause of particle 
acceleration in extragalactic jets [ll| , as a source of high 
energy emission 0: as an explanation of the rapid vari¬ 
ability observed in active galactic nuclei [T^ and many 
others. Therefore, studying magnetic reconnection is of 
great importance, especially considering the fact that the 
relativistic theory of magnetic reconnection is not yet well 
established and its mechanism is poorly understood 0 - 
It should be mentioned that, numerical results of ideal 
relativistic MHD models sometimes show magnetic re¬ 
connection, which is non-physical, since it is caused by 
numerical resistivity and hence depends on the details of 
the numerical scheme and resolution [l^ . 

Therefore, there is a strong interest in developing nu¬ 
merical models for resistive relativistic MHD. However, 
the corresponding governing equations turn out to be nu¬ 
merically very challenging, since the source terms in the 
equations become stiff, especially when the conductiv¬ 
ity is not small [l^. This is the natural consequence of 
the fact that the time-scale of the diffusive effects and the 
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overall dynamical time-scale are of the same order. Thus, 
it is not surprising that the first numerical models for re¬ 
sistive relativistic MHD appeared only in 2006 M and 
2007 [l^. In the latter, the fluxes are computed by using 
the Harten-Lax-van-Leer (HLL) approximate Riemann 
solver and Strang’s splitting technique is used for the 
stiff source terms. Later on, a numerical method which 
uses an implicit-explicit (IMEX) Runge-Kutta method to 
solve the stiff source terms in the equations is proposed 
in Ref.[l3- Also, a unified framework for the construc¬ 
tion of one-step finite volume and discontinuous Galerkin 
schemes for the resistive relativistic MHD is introduced 
in Ref. [l^ . More recently, a different approach has been 
suggested in Ref. 0 , where the method of characteristics 
is used to solve the Maxwell equations. Additionally the 
role of the equation of state in the resistive relativistic 
MHD is investigated in Ref. [2^. 

All the above mentioned models are based on solv¬ 
ing the macroscopic governing equations of the resistive 
relativistic MHD. However, in the last few years, new 
approaches based on lattice Boltzmann (LB) methods 
have been developed to study relativistic hydrodynamics 
(2ll - l^ . Like the other LB methods, they are based on 
a minimal lattice version of the Boltzmann kinetic equa¬ 
tion, where representative particles stream and collide on 
the nodes of a regular lattice with sufficient symmetry to 
reproduce the correct macroscopic equations. The advan¬ 
tages of these models compared to the conventional meth¬ 
ods are their mathematical simplicity, computational ef¬ 
ficiency on parallel computers, and easy handling of com¬ 
plex geometries. 

In this paper, we develop a LB model for relativis¬ 
tic MHD. The model is proposed for resistive MHD but, 
as we will show later, it is robust enough in the ideal 
MHD limit as well. The hydrodynamic part is based on 
the model proposed in Ref. . with several extensions, 
namely to include the contribution of electromagnetic 
fields in the energy-momentum tensor (corresponding to 
adding the Lorentz force and Joule heating in the macro¬ 
scopic equations) as well as to deal with a more gen¬ 
eral equation of state. For the electromagnetic part, i.e., 
solving the Maxwell equations, the LB model for electro¬ 
dynamics proposed in Ref. [2^ is modified and extended 
for coupling with the fluid equations and to include the 
relativistic Ohm’s law. It should be mentioned that in 
the non-relativistic context, there are several LB mod¬ 
els for resistive MHD I25l 4^ . especially for simulating 
magnetic reconnection [28l . l29l| . Our goal is to bring the 
well-known advantages of the lattice Boltzmann schemes 
to the context of resistive relativistic MHD. 

The model is validated using numerical tests for the 
ideal MHD and resistive MHD regimes. In particular, 
the propagation of Alfven waves in high conductivity me¬ 
dia, and the evolution of current sheets in resistive media 
are validated against analytical solutions. Moreover, as 
an application for the model, the magnetic reconnection 
process driven by Kelvin-Helmholtz (KH) instability is 
studied in detail. The KH instability is one of the fun¬ 


damental hydrodynamic instabilities which occurs during 
the shear flow of a uniform fluid, or two fluids with differ¬ 
ent densities. It was discovered independently by Kelvin 
and Helmholtz in the 19th century [^[^. It is believed 
that the KH instability appears in the solar-wind interac¬ 
tion with the Earth’s magnetosphere which can influence 
the magnetic reconnection process that takes place there 
[ 3 ^ . Moreover, the KH instability has been widely inves¬ 
tigated for astrophysical applications e.g., astro phy sical 
jet morphology motion of interstellar clouds and 
clumping in supernova remnants [s^, where in many of 
these phenomena, relativistic and magnetic field effects 
cannot be ignored. Here, we study the KH instability 
as a potential driver of magnetic reconnection. In the 
non-relativistic context this has been discussed in Refs. 
[ 3 ^ 1^ . Here we focus on this phenomenon in the rel¬ 
ativistic context and we are interested in the effects of 
the hydrodynamics parameters, i.e., shear velocity and 
density ratio, as well as the effects of the conductivity on 
the magnetic reconnection rate. 

Furthermore, the results of a three-dimensional simula¬ 
tion of the magnetic reconnection in a stellar flare driven 
by a shear velocity on its photosphere are presented. It 
has been suggested that solar flares are good prototypes 
for stellar flares in relativistic stars like neutron stars . 
Therefore, a solar type initial condition, consisting of a 
potential quadrupole background field and a flux rope 
[3^ . is chosen to study the stellar flare. We show that 
the shear velocity on the photosphere of the star can 
cause the magnetic reconnection to take place between 
the flux rope and the background magnetic field lines. 

The paper is organised as follows: in SecHIl the basic 
equations for resistive relativistic MHD are presented; in 
Sec IHIl the development of a lattice Boltzmann model for 
solving the governing equations is elaborated; in Sec lIVl 
validation tests and the aforementioned applications of 
the model are presented; and finally in SeclVl as a con¬ 
clusion, an overall discussion of the model and the results 
are provided. 

II. THE RESISTIVE RELATIVISTIC MHD 
EQUATIONS 

The equations of motion for resistive relativistic MHD 
can be written in the covariant form as 

d^N>^ = 0 , ( 1 ) 

where = nU^ is the density current with n the 
mass density, {U^) = {c,u)j{u) the four-velocity, u the 
three-dimensional velocity, c the speed of light, 'y(u) = 
i/yr — u^/c^ the Lorentz’s factor, and 

= 0 , ( 2 ) 

where is the total energy-momentum tensor defined 
as the sum of fluid energy-momentum tensor and the 
contribution of the electromagnetic fields, i.e., 

rjifll/ _ rpflV I rpIlV /q\ 

^ ~ Fluid EM’ 
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with 


'^Fluid — i^+P) 




PT 


where Eq. (fT^ is the Maxwell-Ampere and Eq. (US is 
the Maxwell-Faraday equation. The equation for conser- 
(4) vation of current 


where e is the energy density (including rest mass en¬ 
ergy), p is the hydrostatic pressure and tt^'' is the dissi¬ 
pation tensor. On the other hand 


^EM = + \FP^Fp,r^n. 


(5) 


where is the Maxwell electromagnetic tensor defined 
as 




dtPc -l- V • J — 0, 


(14) 


can be obtained by taking the divergence of Eq. (IT^ by 
considering Eg. (ITU)) and Ea. dTT]) as constraints. 

Furthermore, the coupling between the fluid equations 
and Maxwell equations is expressed by Ohm’s law. In 
general, the explicit form of the current four-vector 
depends on the properties of the electromagnetic fields 
as well as the fluid variables. Here we use Ohm’s law for 
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a resistive isotropic plasma as 
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0 

-cB^ 

cBV 

(6) 
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E + u XB — 


{E ■ u)u 


PcU, 


(15) 


where E* and are the electric and magnetic fields, re¬ 
spectively, and eo is the permittivity of free space which 
relates to the permeability of the free space, poy through 
the relation c^eoPo = 1- Note that, throughout this pa¬ 
per, Latin superscripts (subscripts) run over the spatial 
coordinates, while Greek superscripts (subscripts) run 
over the four-dimensional (4D) space-time coordinates. 

In addition to the mentioned hydrodynamics conserva¬ 
tion equations, the governing equations for the electro¬ 
magnetic fields, i.e.. Maxwell equations, also need to be 
considered, which in the covariant form read as 


with a being the conductivity of the plasma. It is worth 
mentioning that in the fluid rest frame Ohm’s law be¬ 
comes 

J = aE, (16) 

and in the limit of cr —> oo one can obtain the well-known 
result for ideal MHD 

E = -UxB. (17) 


= -/roc/'®. 

(7) 

= 0, 

(8) 


where (/^) = (cpc, J) is the four-vector of electric current 
with Pc the charge density and J the three-dimensional 
electrical current while E*^'^ is the Faraday tensor de¬ 
fined as 

E*^'' = (9) 

with the Levi-Civita tensor. 

By choosing an appropriate decomposition of the 
Maxwell tensor, one can show that Eqs. © and © yield 
the familiar Maxwell equations [s^ 


The major difference between the numerical models for 
ideal and resistive MHD originates from the fact that, 
in the ideal case, one can substitute the electric field E 
in all the equations, using a simple algebraic relation, 
i.e., Eg. lfTTl) . and thus one can define the electromagnetic 
induction four-vector [i^. This leads to a con¬ 

siderably simpler and less expensive numerical algorithm, 
compared with the resistive MHD. 

To summarize the governing equations, and by substi¬ 
tuting Ea. (IT5|) into Eas. dTUl) and da, we have 12 equa¬ 
tions, i.e., Eqs. (dD, (©, ((n|), (dSI) and da and 13 un¬ 
knowns, i.e., u, B, E, e, p, n and pc- This system of 
equations will be complete by including the equation of 
state. Here, we consider the ideal gas equation of state 

M 

p = (T-l)(e-nc^), ( 18 ) 


V ■ E = —pc, 
eo 

V • B = 0, 


(10) where, T is the adiabatic index which for ultrarelativistic 
temperatures takes the value 4/3, and for non-relativistic 
temperatures has the value 5/3. 

( 11 ) 



B — V X B = —poJ, 


( 12 ) 


III. LATTICE BOLTZMANN MODEL FOR 
RESISTIVE RELATIVISTIC MHD 


dtB + V X E = 0, 


In this section we describe our lattice Boltzmann model 
(13) to solve the aforementioned governing equations. 
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Relativistic fluid equations 

We start with the description of our LB model to solve 
the equations of motion of the fluid, i.e., Eqs. m and 
©• The relativistic Boltzmann equation, based on the 
Anderson-Witting collision operator is 

P^d^f = -^{f-n), (19) 

TC‘‘ 

with the Maxwell-Jiittner equilibrium distribution func¬ 
tion as 


r'i = Aexp(-p^C/7fcBT), (20) 

where (p^) = {E/cyp) the four-momentum, with the en¬ 
ergy of the particles E defined by 


E = cp° 


mc^ 

V'l - mVc2 ’ 


( 21 ) 


where m is the rest mass and is the magnitude of 
the three-dimensional velocity. Here, / = f{x,p,t) is the 
probability distribution function of particles at time t in 
location x and momentum p, which can be used to calcu¬ 
late both the density current and total energy momentum 
tensor (see below). Note that in LB models, which are 
based on a lattice version of the Boltzmann equation, the 
aforementioned particles are collective degrees of freedom 
(quasi-particles) which can eventually propagate faster 
than light, as long the physical flow speed remains sub¬ 
luminal. Also, r is the single relaxation time, A is a 
normalization constant, ks is the Boltzmann constant 
and T is the temperature. 

As mentioned, our description to solve the fluid equa¬ 
tions is based on the model recently proposed for rel¬ 
ativistic LB [^ . Several modification and extensions 
are implemented in order to firstly, use ideal gas equa¬ 
tion of state (by using Anderson-Witting as the collision 
operator and extending the distribution function) and 
secondly, to include the contribution of electromagnetic 
fields in the total energy-momentum tensor (by adding 
corresponding terms to the distribution function). Fol¬ 
lowing Ref. |23l| . we use changes of variables which read 
as: 




p^/m 


[/M 

5 

Cs 


( 22 ) 


V = cjcs- (23) 

In order to discretize the Boltzmann equation, Eq. (|T^ , 
we write (^“) = {ct/co,Ca), where ct, co and Ca are con¬ 
stants related to the size of the lattice. We use a lattice 
configuration D3Q19 (19 discrete vectors in 3 spatial di¬ 
mensions), which can be expressed as 

r (0,0,0) * = 0; 

Ca = I Co(±1,0,0)fs 1 < * < 6; (24) 

I Ca(±l,±l,0)j^s 7<f < 18, 



where the subscript FS denotes a fully symmetric set of 
points (see FiglT|). 

Using the quadrature rule and knowing that Wi = 1 
the discretized weights can be calculated as (see Ref. [2^ 
for the derivation) 


Wo = 1 + 


4cj 

361c2 


Cj 

„2p2 ’ 

CaCQ 


(25) 


2166c§c2 


(361 - 8c>2) 


(26) 


for 1 < z < 6, and 


Wi = 


1083cg’ 


(27) 


for 7 < i < 18. 

The aforementioned constants can be calculated as 
Ca = — , ct/co = —, co = |(9-2V3). (28) 


In order to increase the numerical stability of the LB 
model (particularly for higher velocities) a flux limiter 
scheme (min mod scheme) is used to discretize the spa¬ 
tial derivative in the streaming term in the Boltzmann 
equation, i.e. pfdafi- The min mod scheme efficiently 
reduces the instability, especially in the presence of dis¬ 
continuities or large gradients. We have 

= rF^-T [htix + Sxea) - h“(x)], (29) 

\dxea\ 


where Bq is a unit vector in the direction of the corre¬ 
sponding spatial coordinate. For the definition of hf{x) 
and more details see Ref. [2^. 

Therefore, the discretized form of the relativistic Boltz¬ 
mann equation takes the following expression: 


fi{x,t + St) - Mx,t) + [hi{x + SxCa) - h“(f)] = 
Ct dx 

TCt 

Ct ^ 

(30) 

where the last term is the bulk viscosity term with 


fO i = 0] 
1 aSx i ^ 0, 


(31) 


where a is a small constant. The exact value of a for each 
simulation is reported in Sec lIVI In fact, this small value 
of bulk viscosity is sufficient to stabilize the numerical 
procedure. A central finite difference scheme is used to 
calculate the second order derivative. 

One can notice that to solve the discretized Boltzmann 
equation. Eg. (15(11) . the discretized equilibrium distribu¬ 
tion function is required. The discretized equilibrium 
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distribution function to recover has been proposed 

in Ref. [ 2 ^. To include the electromagnetic contribution 
to the energy-momentum tensor, i.e., additional 

terms need to be added to the discretized distribution 
function. Let us elaborate the contribution of electro¬ 
magnetic fields in the energy-momentum tensor by pro¬ 
viding the components of using Eqs. ([5]) and ([6|). 
Thus, we have 


^00 _ 

-^EM — 




c^B^) 


(32) 


rpOi _ 

^EM — 


— {E X B)\ 

^J.qC 


(33) 


-E^E^ - c^B^Bi + -{E‘^ + , 

(34) 

where E'^ and B^ are the magnitude of the electric and 
magnetic fields, respectively, and is the Kronecker 
delta. Note that, adding these terms to the total energy- 
momentum tensor corresponds to adding the Lorentz 
force and Joule heating to the macroscopic equations. 
The following expression shows the complete discretized 
equilibrium distribution function which recovers the to¬ 
tal energy-momentum tensor with the ideal gas equation 
of state 


'EJm — ^0 


/r = 7(e+p)3«^. 1 + 


3(r — l)(e — nc^) — e 
(r — l)(e — nc^) + e 


361 [e — (T — l)(e — nc^)] 
33[(r-l)(e-nc2)-e] 




+ c^clx^X^ + clclx^r + (|^ - • X) 


VCq 


+ +{cinx^^+{cinxr 

- ^(x • x)] I + + E^) 

+ i[c\B-B) + E-E]-^[{BxE)-Ca] 

+ y [c^{B ■ caf + {E • ca)2] - ^ [clclE^Ey 
+ clclE^E^ + clclE^E^) + c\clclB^By 


EclclB^B^ ^clclB^B^)^, 

(35) 

where the second curly bracket is the contribution of the 
electromagnetic fields. Note that the second and third 
term in the first curly bracket are the contribution of 
the ideal gas equation of state which goes to zero in the 
ultrarelativistic limit (T = 4/3 and e ^ nc^) where the 
equation of state becomes e = 3p. 

As mentioned. Eg. (1501) can be used to solve the equa¬ 
tions for the conservation of the total energy-momentum 



FIG. 1. The D3Q19 lattice configuration for the LB model to 
solve the relativistic fluid equations. 


tensor, Eq.(l2]). To solve the equation for the conserva¬ 
tion of number of particles, Eq. CD, a separate distribu¬ 
tion function based on the model proposed in Ref.ji^ 
is considered. Therefore, we add an extra distribution 
function, gi, which follows the dynamics of the Boltz¬ 
mann equation given by Eq. (1301) . without the \i coeffi¬ 
cient term and by substituting (C/^X^/^^) by unity. The 
corresponding modified equilibrium distribution function 
for the cell configuration given in Ea. (l24)) is given by: 


97 = w[ni{u) ( y -b 3(cY ■u) + -{ca- uf 


with 


, _ 1 , _ 3 1 

“"o 10’ 10 QcV 


for 1 < f < 6, 


1 3 


Wi = 


I2cl 40’ 



(37) 


(38) 


for 7 < i < 18, where w[ are the respective discrete 
weights. Note that computing the macroscopic variables 
for the fluid, is not any more straightforward. We shall 
elaborate this issue later on. 


Maxwell equations 

Now that the LB model for solving the fluid equations 
is discussed, let us explain our LB model for solving the 
governing equations of electromagnetic fields, i.e., Eqs. 
C2D,C3D, and d, with (TTKll as Ohm’s law. Our scheme 
is based on a 3D LB model for solving the Maxwell equa¬ 
tions proposed in Ref. [ 2 ^, where several modifications 
are required to couple it to our solver of the fluid equa¬ 
tions (mainly by modifying the distribution functions) 
as well as to use it for relativistic MHD (by using the 
relativistic Ohm’s law). 
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For this purpose, we use a cubic regular grid with 13 
velocity vectors (D3Q13), where four auxiliary vectors 
are assigned to each of the vectors (two for the electric 
field and two for the magnetic field) for calculating the 
magnetic and electric fields. Note that, since Maxwell 
equations are first order, non self-interactive differential 
equations, the symmetry requirements are less strict than 
the symmetry requirement for the fluid dynamics equa¬ 
tions. Therefore, compared to the lattices for fluid dy¬ 
namics, less velocity vectors, e.g., 13 velocity vectors, 
are required. Indeed, one could also use the D3Q19 for 
the Maxwell equations, but it would unnecessarily spend 
more computational resources. 

A simple streaming-collision evolution for the distribu¬ 
tion function is considered as 

'^h 

(39) 

and 

hg (x, t-I-(5t) - /ig(x, t) = - — [hg(f, t) - /ip(x, t)], (40) 

Th 

where (x, t) and (x, t) are the equilibrium dis¬ 
tributions to be defined later. Here i = 0,1, 2, 3 indicates 
the direction of the vectors, p = 0,1, 2 shows the plane 
where the vectors lie, and j = 0,1 shows each of the two 
auxiliary vectors for the electric or magnetic field. Thus, 
there are four directions on three planes which gives 12 
vectors, and including the rest vector, in total we have 
13 vectors. These vectors (except the rest vector) lie on 
the diagonals of each plane and point to the edge-centers 
of a cube, so we can write the components as 

■u® = 2c{cos {2i + l)7r/4, sin {2i + l)7r/4,0} , (41) 

vl = 2c{cos {2i + l)7r/4,0, sin {2i + l)7r/4} , (42) 

cf = 2c (0, cos {2i + l)7r/4, sin {2i + l)7r/4} , (43) 

in addition to the rest vector, i.e, vq = (0,0, 0). 

The distribution functions propagate with these vec¬ 
tors from cell to cell. Note that, unlike the LB models for 
fluid dynamics, these vectors do not represent the veloc¬ 
ity of any particle. Associated to each velocity vector cf 
there are two electric auxiliary vectors eff and two mag¬ 
netic auxiliary vectors 6A, which are used to compute the 
electromagnetic fields. These vectors are perpendicular 
to uf’. However, lies on the same plane as cf, while 

&A lies perpendicular to this plane. More accurately, we 
define them as (see Figl2|) 

sP - ^ (44) 

and 


where (i)mod4 is a function that gives the remainder on 
the division of i by 4. To these vectors we shall add the 
null vectors, i.e., eb = (0,0,0) and bo = (0,0,0). This 
means that there are 13 different electric vectors and 7 
different magnetic vectors. 

In order to solve the Maxwell equations by the LB 
model we can write Ampere’s law (Faraday’s law) as time 
derivative of electric (magnetic) field plus the divergence 
of an antisymmetric tensor [^. We also consider the 
term(— /igj) in Ampere’s law (right hand side of Ea. dT^ l 
as an external force. Therefore, the macroscopic fields 
can be computed as 

(«) 

p—0 j—0 

s = EEE'‘t%. («) 

i—0 p—0 j—0 

and 

3 2 1 

Pc = ho + ^^^hPy (48) 

z^O p—0 j—0 

Note that the effect of the external force still needs to be 
considered to get the correct electric field and E' is the 
electric field before considering the external force. 

It can be shown that to recover the Maxwell equations, 
for the current model, ^ should be considered. Un¬ 
like the LB models for fluid dynamics, this value for the 
relaxation time does not lead to numerical instabilities, 
because of the linear nature of the Maxwell equations. 
For the case oirh = \ the external force in Ampere’s law 
can be included in a rather simple way, and E becomes 
(see Ref. [^ 1 

E = &- jpoc^l (49) 

where, according to the Ohm’s law, Ea. lfTKll . J is a func¬ 
tion of E. By substituting Ohm’s law in Eg. lHHl) we ob¬ 
tain a system of three equations and three unknowns {E^, 
E^ and E^), which can be solved analytically. Having the 
values of each component of the electric field, J can be 
calculated using Ohm’s law, consequently. More discus¬ 
sion about computing the macroscopic variables shall be 
provided later. 

The discretized equilibrium distribution function to re¬ 
cover the correct Maxwell equations reads as follows 

hg®'’jx, t) = • J -f -^E ■ -k ■ 1^^-, (50) 


I 

2^)2' 


X et„ 


(45) 


and 


= Pc. 


(51) 


7 



FIG. 2. The configuration of the auxiliary vectors for the LB 
model to solve the Maxwell equations. 


A Chapman-Enskog expansion shows that the current 
model recovers Ampere’s law, Faraday’s law and current 
conservation, Eq. m- The latter follows from the evolu¬ 
tion of ho- 

The divergence free condition for the magnetic field, 
Eq. (HID , can be treated as a constraint on the initial con¬ 
dition, since by taking the divergence of Faraday’s law 
one can show that the time derivative of V • i? is always 
zero. Therefore, if V • B = 0 is set for the initial condition 
it will hold for later times as well. The same is true for 
the Gauss law, Ea. dTOl) . using Eqs. (HU and ([Mil . 


Coupling between fluid and electromagnetic fields 


Having explained the appropriate solvers for fluid 
equations and Maxwell equations, we next discuss how to 
compute the macroscopic variables. As mentioned, the 
model of Anderson-Witting is used for the collision term 
in the solver for the equation of conservation of energy- 
momentum. Hence, according to the Landau-Lifshitz de¬ 
composition [ 3 ^ the macroscopic variables can be cal¬ 
culated by solving an eigenvalue problem resulting from 
multiplying the relation for the energy-momentum tensor 
by the covariant four-vector velocity. Using the definition 
of the total energy-momentum tensor, Eqs. ([3]),(HI) and ([S]) 
with the help of the relation = 2{c^B^ — E'^) we 

get 





E;2) U'', (52) 


Here, [e -I-— E^)] and are the largest 
eigenvalue and corresponding eigenvector of the ten¬ 
sor [T^ — eoF^^F/yp], respectively. The total energy - 
momentum tensor can be calculated using the relation 
On the other hand, the tensor 
F^^Fyp depends on E and B. As mentioned, for cal¬ 
culating E the value of the external force, which de¬ 
pends on the velocity, is required. However, the value 
of the velocity is not yet computed. To solve this 
problem, we use the value of the electromagnetic fields 
from the previous time step to calculate the tensor 
\Tjf — e^F^PEvp]. The largest eigenvalue and correspond¬ 
ing eigenvector (velocity) can be calculated numerically 


using the power method. Knowing the velocity, one can 
calculate the density using the first order moment rela¬ 
tion, i.e., = nU^ = After that, B is 

calculated using Eg. dTTl) . where, the are obtained 
using the LB equation for the Maxwell equation, i.e, 
Eq. (I39|) . in which the (Eq. (I50|) l are known from 

the previous time step. Having the velocity and magnetic 
field, E can be computed using Eo. ddfip and including 
the external force as described before. Thus, it is easy 
to compute J using Ohm’s law. Considering the eigen¬ 
value [e + ^ (c^B^ — E^)], it is possible to compute e and 
through the equation of state Eq. (ITS)) one can compute p. 
Finally, can be calculated directly from Eg. (ITS)) . All 
the 13 unknown variables for each cell can be computed 
in this way. Note that using the values of electromagnetic 
fields from the previous time step to calculate the ten¬ 
sor, leads to an error which goes to zero as the time step 
(St) decreases. In fact, in the next section our numerical 
results show that the error is ignorable. 

Finally, since two separate solvers for fluid equations 
and Maxwell equations are used in this model, we need to 
make sure that the time evolution of the two solvers are 
consistent with each other, i.e. the distribution functions 
for the fluid and the electromagnetic solvers must evolve 
simultaneously in time. For the electromagnetic solver 
the time evolution, i.e. St depends on the value of the 
spatial spacing, i.e., Sx through the relation Sx/St = v^c, 
whereas for the fluid solver Sx/dt can be chosen freely, 
as long as the stability of the numerical model is not 
impaired. Since the value of Sx is the same for both 
solvers, to ensure the same value of St for both solvers, 
the numerical value of the velocity of light is adjusted 
accordingly, and used for both solvers. Also, note that, 
in order to describe signals that move at the speed of 
light, e.g. electromagnetic waves, the numerical quasi¬ 
particles in the LB model should move faster than light 
(according to the relation Sx/St = v^c). Nevertheless, 
in the continuum limit the differential equations that our 
model reproduces, which characterize the physical sys¬ 
tem, never violate the principle of relativity. 


IV. TEST SIMULATIONS AND APPLICATIONS 

In this section, we present some numerical tests in or¬ 
der to validate our numerical model along with some ap¬ 
plications for the resistive relativistic MHD. More specifi¬ 
cally, test simulations for the propagation of Alfven waves 
and the evolution of a self-similar current sheet are con¬ 
sidered, and as applications for the resistive relativistic 
MHD LB model, we study the magnetic reconnection 
driven by the Kelvin-Helmholtz instability and present 
the results of a 3D simulation of magnetic reconnection 
in a stellar flare due to the shear velocity in the photo¬ 
sphere. The purpose of the first test simulation (Alfven 
wave), is to validate the numerical method in the limit of 
ideal MHD, while the second test is to validate the model 
recovering the correct dynamics in the resistive regime. 









X 


(b) 

FIG. 3. (Color Online) Results of the numerical test for the 
propagation of Alfven waves in the limit of ideal MHD (cr = 
10®) for (a) magnetic field in 2 -direction and (b) electric field 
in x-direction. For both plots, the numerical results are shown 
at t = 1 (blue circle symbols) and at t = 1.5 (red square 
symbols). The dashed black line shows the initial condition 
and solid black lines are the exact analytical solutions at the 
corresponding time. 


In the following simulations numerical units are used and 
/io = 1 is considered. 


region of the initial wave it is assumed that ft = 0 and 
= 0. Inside the region of initial wave we have 

= tiaBq sin r27r(3a;^ — 2x^)1 , x* = — —(53) 

Xi- Xq 

and 


= uy = 0, 

Bo 

where the Alfven velocity is defined as 


(54) 


= 


2Bl^ 


e + p + Bl{l + r]\) 


1 + 1 1 - 


2r^ABl 


+ P + Bl{l + ri\) 


(55) 


Here we use the value of rjA = 0.118591 which gives va = 
0.40785. Ohm’s law for ideal MHD, i.e. Eg. lflTl) is used to 
calculate the initial electric field and Ea. dTUl) to measure 
the initial value of pc- To achieve the ideal regime a very 
high conductivity (cr = 10®) is considered. The equation 
of state, Eq. dUD, with r = 4/3 is used. The domain is 
discretized using 400 cells and we set 5x/5t = v^, which 
gives c = 1. The value of r for the fluid LB model is set 
to 1 and a = 0.1. Open boundary conditions for each 
cell at the left and right boundaries are implemented by 
copying the distribution functions from the neighboring 
cell in the direction perpendicular to the boundary. 

The simulation runs until t = 1.5 and the results are 
presented in Eigl^l In the limit of ideal MHD the gener¬ 
ated wave should travel with the Alfven velocity without 
any distortion. Here we have compared our numerical re¬ 
sults with the exact analytical solution of the ideal MHD 
at two different times, i.e., t = 1.0 and t = 1.5. The 
results are presented for the magnetic field in 2 -direction 
(Fig 3(a) I and the electric field in ^-direction (Fig 3(b) |. 
Figli shows that at each of the considered times and for 
both B^ and E^, very good agreement is observed be¬ 
tween the numerical and analytical results. This simu¬ 
lation takes ~ 320 ms on a single core of an Intel CPU 
with 2.40 GHz clock speed. Note that, the analytical 
results are obtained by simply shifting the initial wave 
by the Alfven velocity. Apart from validating the nu¬ 
merical method, this test shows the ability of the model 
to deal with high conductivity (low resistivity) regimes, 
recovering the ideal MHD limit. 


Propagation of Alfven wave 

This test deals with the propagation of an Alfven wave 
along a uniform background field in the limit of ideal 
MHD. The initial condition is the same as in Ref. [l5l| . 
We set n = 1.0, p = 1.0, B^ = Bq = 1.0 and By = 0.1 
and we consider a one dimensional domain defined in the 
range of — 1 < a; < 1, where the initial wave is located at 
xo < X < Xi, with xo = —0.8 and Xi = 0. Outside the 


Evolution of self-similar current sheets 

After validating the model for the ideal MHD case, we 
consider here a test problem in the resistive case for which 
the evolution of a current sheet is investigated. We as¬ 
sume that the magnetic pressure (H^/2) is much smaller 
than the plasma pressure (p), so that the fluid is not af¬ 
fected by the evolution of the current sheet and changes 
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in the magnetic field. We know that when the magnetic 
field changes its sign within a thin layer a current sheet 
forms. Thus, for our case, we assume that the magnetic 
field has only a tangential component B = (0,i3^,0), 
where = B{x,t) changes sign within a thin current 
sheet of width Al. If the fluid is set initially to equilib¬ 
rium, by considering a constant pressure in the domain, 
the evolution of the current sheet becomes a diffusion pro¬ 
cess. By assuming that the diffusion time-scale is much 
longer than the light propagating time-scale, we can ne¬ 
glect the displacement current {dtE) in Ampere’s law. 
In the rest frame, by inserting the relation J = aE in 
Ampere’s law, using Faraday’s law, and plugging in the 
mentioned one-dimensional magnetic field of the current 
sheet, one gets 

dtBy - -dlBy = 0. (56) 

a 

As the diffusion process continues and the width of the 
current sheet becomes much larger than the initial width 
(A/), the expansion becomes self-similar and the analyt¬ 
ical solution has the form 

B{x, t) = Boeii . (57) 

where Bq is the magnetic field outside of the current sheet 
and erf is the error function. The above equation de¬ 
scribes the evolution of a current sheet, providing dtE is 
ignorable. 

For the numerical test a domain of—1.5 < x < 1.5 is 
discretized using 100 cells, where open boundary condi¬ 
tions are considered for the left and right boundaries. 
The initial values p = 50,n=l,zi = £' = 0 and 
B^ = = 0, are considered, and the initial B^ is com¬ 

puted using Eg. (1571) at t = 1 with Bq = 1. In the nu¬ 
merical model F = 4/3, 5x/5t = v^, t = 1 and a = 0.1 
are considered. The simulation runs until t = 8 and the 
results are compared to the analytical results at t = 9 
(since the initial condition is assumed to be at t = 1). 
Two values of uniform conductivity, a = 100 and a = 50, 
are considered and the results of the comparison with the 
analytical solution for both cases along with the initial 
conditions are presented in FigU] One can see that the 
numerical and analytical results are almost indistinguish¬ 
able and that the current sheet in a domain with cr = 50 
diffuses faster than the domain with tr = 100 due to the 
higher resistivity. This test validates the resistive part 
of the numerical model and shows the capability of the 
model for simulating resistive problems far from the ideal 
MHD limit. The above simulation, for cr = 100, takes 
^ 70 ms on a single core of an Intel CPU with 2.40 GHz 
clock speed. 

To check the convergence of the model, we implement 
the same current sheet simulation with ct = 100 for dif¬ 
ferent grid resolutions. Fig IS] reports the error versus the 
number of cells in a log-log plot, and the slope of the fit¬ 
ted line (blue solid line) shows that the model is second 



FIG. 4. (Color Online) Results of the numerical test for the 
evolution of a current sheet in the resistive regime. Dashed 
green line (upper curve at a; > 0) and dashed red line (sec¬ 
ond upper curve at a: > 0) show the initial condition of the 
current sheet for cr = 100 and a = 50, respectively, which 
corresponds to the analytical solution. Eg. (1571) . at t = 1. The 
green triangle and red square symbols show the result of the 
numerical solution at t = 9, for a — 100 and cr = 50, respec¬ 
tively. The solid black lines show the exact analytical solution 
at this time. 


order as we expect for a lattice Boltzmann scheme. The 
error is calculated as follows: 

/ N \ 1/2 

U E(^num - > (58) 

where N is the number of cells, and B^^^ and 
are the numerical results for the cases with N and 200 
grid resolution, respectively. Due to the fact that the an¬ 
alytical expression, i.e., Eq. (I57|), is only approximately 
accurate, a fine mesh (200 cells) is used as a reference to 
calculate the error. In other words, for each grid resolu¬ 
tion, the numerical result at each position is compared to 
the numerical results of the fine resolution. Additionally, 
as explained in the introduction, one of the advantages 
of using lattice Boltzmann methods is its simplicity and 
efficiency on parallel computers. To show that this also 
holds for our model, we use a simple openMP paralleliz¬ 
ing method and simulate the 3D extension of the afore¬ 
mentioned current sheet problem with 150 x 150 x 150 
cells and cr = 100, until 5 time steps. The preliminary 
resulting speedup for a few number of threads is shown 
in the inset of EiglS] where one can see that even for a 
straight forward parallelizing method, a satisfactory level 
of efficiency is achieved. However, for a thorough study 
of the parallelizing efficiency, a much larger number of 
threads and more advanced parallelizing methods should 
be experimented, which is an interesting topic for future 
research. 
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FIG. 5. Results of the convergence test for the simulation of 
the evolution of a current sheet in the resistive regime with 
a = 100 in a log-log plot. Symbols show the numerical error 
and the solid blue line is the fitted line with slope —1.9668, 
which shows a second order convergence. In the inset the 
speedup obtained by using the openMP parallelizing method 
is reported versus the number of threads. The dashed line 
shows the ideal speedup. 



Magnetic reconnection driven by Kelvin-Helmholtz 
instability 


After validating the numerical model, we now study 
the magnetic reconnection process driven by the Kelvin- 
Helmholtz (KH) instability. As mentioned before, the 
magnetic reconnection is a process where magnetic lines 
change their topology. At the place where the magnetic 
lines reconnect, usually a null point forms where the mag¬ 
netic field vanishes. There are several theoretical descrip¬ 
tions for the magnetic reconnection inclu ding the well 
known Sweet-Parker [13, Petschek |46| models. 

The Sweet-Parker model is based on the discussion of 
pressure balance in the reconnection region, where the 
reconnection region is assumed to be dominated by dif¬ 
fusion while the outside region is assumed to be ideal. 
Due to the magnetic diffusion, the plasma is driven into 
the current sheet (inflow) with velocity Uin in the direc¬ 
tion perpendicular to the length of the current sheet. 
The conversion of magnetic energy within the current 
sheet thrusts the plasma out with the velocity Uout in 
the direction of the length of the current sheet. It is 
shown that the reconnection rate R, which is defined as 
the ratio between Vm and Uout is proportional to S~'^, 
where S = ^qLvao' is the Lundquist number (magnetic 
Reynolds number), with L as characteristic length. The 
reconnection rate in this model is usually small due to 
the high aspect ratio of the reconnection region (which is 
proportional to the inflow and outflow velocities assum¬ 
ing the incompressibility condition). 

In the Petschek model, it is assumed that the magnetic 
energy can be liberated not only in the current sheet but 


0.5 1.0 1.5 

FIG. 6. (Color Online) Snapshots of the density for the KH 
instability at times (a) t = 0.0 (initial condition) (b) t = 3.31 
(c) t = 6.62 (d) t — 15.47, for the case with a = 100, An = 
1.8, Uo = 0.6c. The white lines show the magnetic field lines. 



FIG. 7. (Golor Online) Reconnection rate versus time for the 
case G = 100, An = 1.8 for different values of Uq. In the inset 
the results are shown for two different values of An =1.8 and 
An = 0 for a = 100. 
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also as pairs of slow shocks which stem from the edge 
of the sheet. Therefore, the reconnection region can be 
smaller than the one for the Sweet-Parker model, which 
can lead to higher reconnection rate. In this model R is 
proportional to (10 5)“^. 

As discussed, magnetic reconnection is believed to 
have prominent effects in many high energy astrophysical 
events and therefore it should be strongly influenced by 
relativistic effects. Although the mechanism of relativis¬ 
tic reconnection is not well understood, recent theoretical 
and numerical studies of the relativistic Sweet-Parker and 
Petschek models show the same proportionality relation 
between R and S [11,143. In the numerical simulations, 
it is important how one triggers the reconnection process. 
If a local increase in the conductivity is used to trigger 
the reconnection, Petschek type reconnection is observed 
with “x-type” null point, while when a perturbation in 
the magnetic potential is applied, Sweet-Parker type re¬ 
connection is observed with “y-type” null point [11] . This 
is similar to the non-relativistic numerical results. Here 
instead of using either of the mentioned ways, we use 
a hydrodynamic instability to trigger the reconnection. 
In particular, the KH instability is chosen because of its 
wide range of applications in astrophysical events as dis¬ 
cussed in the introduction section. Therefore, two di¬ 
mensional simulations of the KH instability in a Harris 
like current sheet are performed for different values of 
shear velocity, density ratio and conductivity. A domain 
of —0.5 < a; < 0.5 by —0.5 < ?/ < 0.5 is discretized us¬ 
ing 512 X 512 cells. The following initial conditions are 
considered: 

Mx =tanh , n = no-ftanh [-V (59) 

2 Va/ 2 Va/ 

where I/q defines the shear velocity, An = (uup — Udown) 
is the density difference between the upper {y > 0) and 
lower {y < 0) part of the domain, uq = 1 , p = 1 and 
a = is considered. Furthermore, 

H^ = Hotanh(0, (60) 

while Ry = = 0 and Bq = 0.06 is considered. Note 

that high values of magnetic field will stabilize the KH 
instability due to the magnetic tension of the background 
magnetic held. The KH instability is triggered by a per¬ 
turbation in the velocity in y-direction, namely 

Uy = itpert sin (kx) exp > (^1) 

with Upert = 0.01 as the perturbation amplitude, fc = 27r 
for a single mode perturbation, and b = 10i5a:. For the 
left and right boundaries {x = ±0.5) periodic boundary 
conditions are considered while for the upper and lower 
boundaries (y = ±0.5) open boundary conditions are im¬ 
plemented. The simulation is performed for different val¬ 
ues of Uq = 0.6c, 0.4c, 0.2c and a = 100,80, 60,40, 20 for 
two values of An = 1.8 and An = 0 where the latter 



0.96 1.0 1.04 

FIG. 8. (Color Online) Snapshots of the density for the KH 
instability at times (a) t = 0.0 (initial condition) (b) t = 3.31 
(c) t = 6.62 (d) t = 15.47, for the case with a — 100, An = 0, 
Uo ~ 0.6c. The white lines show the magnetic field lines. 


corresponds to the case with initial uniform density. For 
the numerical simulation F = 4/3, Sx/St = 2.5-\/2, r = 1 
and a = 0.1 are considered. The initial electric held can 
be simply computed using Faraday’s law (dropping the 
time derivative term because of the stationary state) and 
Ohm’s law, knowing the fact that the initial velocity and 
magnetic helds are in the same plane. Additionally, for 
the open boundary condition, and to ensure the diver¬ 
gence free condition of the magnetic held, the normal 
component of the magnetic held is adjusted in order to 
have V • H = 0 at the boundaries. This will slightly im¬ 
prove the numerical stability of the model. Note that 
for the electric held, the same idea (setting the normal 
electric held on the boundaries according to V ■ E = pc) 
does not affect the numerical stability. Therefore, for 
each cell, the open boundary condition for the electric 
held is implemented like for the other variables, i.e., by 
copying the distribution functions from the neighboring 
cell in the direction perpendicular to the boundary. 

The snapshots of the density for the case with tr = 100, 
An = 1.8 and Uq = 0.6c are presented in Figl6]for dif¬ 
ferent time steps. One can see the evolution of the in¬ 
stability in time where after an initial linear growth, a 
non-linear stage takes place which leads to the penetra- 
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FIG. 9. (Color Online) Reconnection rate versus time for the 
case An = 1.8 and Uo = 0.2c for different values of a. 


tion and mixing of the lighter and heavier fluids, where 
the characteristic structure of the KH instability forms. 
The reconnection is triggered by the instability and oc¬ 
curs in the initial current sheet, where the magnetic field 
changes sign. As the instability evolves, the location of 
the reconnection null point changes. Since no external 
force is implemented, the instability finally smoothens 
out due to the dissipation in the system. 

We are interested in the reconnection rates of the con¬ 
sidered cases. Since the initial condition is not in the rest 
frame, it is not practical to find the reconnection rate 
based on the definition mentioned before (ratio between 
inflow and outflow velocities). Instead we use 0 


m 



(62) 


i.e., the generated out-of-plane component of the electric 
field (AA* = Eq being the initial out-of-plane 

electric field) normalized by Bq at the null point, which 
shows the rate at which the magnetic flux is convected to 
this point (see Ref. [H). Here we investigate the effects 
of different parameters on the reconnection rate. First 
we study the effects of the parameters related to the KH 
instability, namely Uq and An. Fig. [7] shows the results 
for the reconnection rate versus time for different values 
of Uq, when cr = 100 and An = 1.8. It is shown that an 
increase in the magnitude of the shear velocity reduces 
the reconnection rate. This is due to the stabilizing ef¬ 
fects of the shearing velocity on the tearing instability 
, an MHD instability that appears in connection with 
sheared magnetic field. Also note that the bumps that 
appear in the reconnection rate at high shear velocities 
is due to the transition of the instability into its station¬ 
ary state, where one can appreciate that, for higher shear 
velocities, this transition occurs at earlier times. 


Additionally, in the inset of Figl?! for the case with 
a = 100, the effects of different values of An is shown. 
One can notice that changing the value of An from 1.8 to 
0 (initially uniform density) has negligible effects on the 
reconnection rate, for different magnitudes of the shear 
velocity. This is despite the fact that the hydrodynamics 
of the system for the initially uniform density is quite dif¬ 
ferent from the case with inhomogeneous densities. FigjS] 
shows the snapshots of the density for the case An = 0, 
a = 100 and Uq = 0.6c, where the well known “cat’s eye” 
structure of KH instability for the case with initially uni¬ 
form density can be recognised. Comparison between 
FiglHand Figl8]shows that, for the case with the initially 
uniform density, the results are symmetric and because 
of the form of the initial perturbation, the location of the 
null point is always at the boundary, unlike the previous 
case, where the location of the null point changes with 
time. 

Another parameter that we are interested in studying 
is the conductivity tr. FiglH] shows the results of the re¬ 
connection rate versus time for the case with Uq = 0.2c 
and An = 1.8 for different values of the conductivity. As 
shown, the reconnection rate increases faster in time and 
reaches a higher value for lower conductivities (higher re¬ 
sistivity). The fact that the reconnection rate increases 
by increasing the resistivity is also expressed in the model 
of Sweet-Parker and Petschek. The interesting point is 
to inspect the exact relation between the reconnection 
rate and the resistivity. As mentioned before, in the 
Sweet-Parker model, R is proportional to a~^ (assum¬ 
ing constant Alfven velocity) and in the Petschek model 
R is proportional to (Incr)“^. To compute this propor¬ 
tionality relation for our results, the reconnection rates 
R{t = ti) at the final time ti = 23.2 are used. The 
results are shown in FigHUl where, as expected, Rq de¬ 
creases with increasing a. The blue dashed line, and 
the red dashed-dotted line in Figlin] are the best fitting 
curves for the data using the proportionality relations 
suggested by Sweet-Parker and Petschek models, respec¬ 
tively. One can see that the results clearly do not follow 
the Petschek scaling law while match the Sweet-Parker 
scaling law very closely. 


Three-dimensional magnetic reconnection in a 
stellar flare 

In this section, we show the results of the 3D numeri¬ 
cal simulation of magnetic reconnection in a stellar flare, 
which is driven by a shear velocity on its photosphere. 
Thus, a domain of —3 < a; < 3 by —3 < ?/ < 3 by 
0 < z < 6 is discretized by 256 x 256 x 256 cells. The 
conhguration of the initial condition is chosen to mimic 
the arcade and the flux rope of a stellar flare [5l[. The 
total potential held (flux function) is defined as 


-ip = tpb + i’i + tpi- 


(63) 
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FIG. 10. (Color Online) R{t) at time t = 23.2 for differ¬ 
ent values of a based on the numerical results (red symbols). 
Best fitting curves using the proportionality relation for the 
Sweet-Parker (dashed-dotted blue) and Petschek (dashed red) 
models are shown. 


This configuration consists of a background magnetic 
field that is produced by four line electrical currents (just 
below the photosphere), which determines ipb and an im¬ 
age electrical current (below the photosphere), which de¬ 
termines tpi. This background magnetic field yields a null 
point above the photosphere. Additionally, a line elec¬ 
trical current contained within the flux rope with finite 
radius, which determines "0;, is added to the null point. 
Note that if we do not consider the flux rope, the mag¬ 
netic field configuration is a potential quadrupole field. 

The potential fields generated by these currents are 
defined as follows: 

V’f) = 

[(x + 0.3)2 ^ ^ 0.3)2][(a: - 0.3)^ + {z + 0.3)^] 

[{x + 1.5)2 + {z + 0.3)2][(x - 1.5)2 + {z + 0.3)2] ’ 

(64) 

^i = -Y^og[x^ + {z + hf], (65) 





r < To, 


^ - To log To + ro log r, r > Tq, 


( 66 ) 


where r = [x'^ + {z—h)'^]^^^ is the distance from the center 
of the flux rope, h is the height of the flux rope, which is 
set to 2.0, ro is the radius of the flux rope, which is set 
to 0.5, and (±1.5,—0.3) and (±0.3,—0.3) are the (x,z) 
positions of the four line currents. Here, Cb represents the 
strength of the background magnetic field, which is set 
to 0.2534 [5l|. The resulting magnetic field is calculated 
by taking the curl of the total potential field, i.e., B* = 
V X ijjey, and in order to adjust the magnitude of the 
magnetic field, each component of the magnetic field is 
multiplied by where B*^^^ is the maximum 



(a) (b) 
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FIG. 11. (Golor Online) Snapshots of the 3D magnetic re¬ 
connection in a stellar flare due to the shear flow at times 
(a) t = 0 (initial condition), (b) t = 19.86 and (c) t = 43.02. 
The colors of the magnetic lines show the magnitude of the 
magnetic field, where blue to red show low to high values. On 
the outer boundary surfaces of the domain the colors indicate 
the value of the vorticity which is indicated in the colorbar. 


of the computed magnetic field in the domain and Bq = 
0.06 is considered. To satisfy the force balance inside 
the flux rope, a magnetic component should be added 
in the direction perpendicular to the magnetic field, i.e., 
y-direction, which has the following form: 


By = 


(Bo/i?^.ax)y2(l-^), r<ro. 


(67) 


r > r-Q. 


One can see the illustration of the initial configuration of 
the magnetic field in Fig. 11(a) The initial velocity is 
set to zero, and the initial density and pressure are con¬ 
sidered to be uniform and equal to 1.0 everywhere in the 
domain. In order to compute the initial values for the 
electric field and electrical current density, the numeri¬ 
cal code is used in an iterative process where the value of 
magnetic field and velocity is fixed at each time-step. The 
converged values of the electric field and electrical current 
provide the initial condition for these variables. Open 
boundary conditions are considered for all the bound¬ 
aries, and as mentioned before, the normal component of 










14 



force, two vortices form, which rotate in the same di¬ 
rection and therefore give rise to a shear velocity in the 
middle of the xy plane, i.e, the photosphere plane. This 
shear velocity finally results in a “cat’s eye” structure 
in the photosphere (see Fig fTTT cll similar to the results 
of the 2D KH instability with uniform initial density in 
FiglU The shear velocity starts to twist the foot of the 
background magnetic lines and later the upper parts of 
the background magnetic lines. As a result at some point, 
the twisted background magnetic lines and the flux rope 
magnetic lines take opposite directions. This is where the 
current sheet forms and the reconnection between these 
two sets of magnetic lines takes place (see Fig fTTT b')'). At 
late times. Fig. HHc), most of the flux rope magnetic 
lines reconnect with the background magnetic lines, and 
only a small part of the flux rope can reach the opposite 
surface. This process can be appreciated more clearly in 
the 2D projections on the xz plane, which are provided 
in Figllll For instance, a starting configuration of the re¬ 
connection can be observed in Fig ll2f bb where one can 
see that the background magnetic lines reconnect to the 
flux rope lines and change their topology. 


V. CONCLUSIONS 


FIG. 12. (Color Online) Projection of the results presented 
in Fig llll onto the xz plane. Times and colors are the same 
as Figim 

the magnetic field is adjusted on each boundary to satisfy 
the divergence free condition of the magnetic field. 

To set-up a shear velocity in the photosphere, two adja¬ 
cent vortices rotating in the same direction are produced 
by applying proper external forces. Therefore, an exter¬ 
nal force is applied to both sub-domains of —3 < a; < 0 
by —3 < ?/ < 3 by 0 < z < 0.2343 as well as 0 < x < 3 by 
<y<3by0<2;< 0.2343, with the following form; 

= Fo sin (7rx*/(3)) cos (7ry/6), (68) 

= Fq cos (7rx*/(3)) sin {iry/d), (69) 

where F^ and F^ are the components of the external 
force in the x and y directions, respectively, Fq is the 
magnitude of the applied force, which is set to 0.05 and 
X* = X for the sub-domain —3 < x < 0 and x* = x — 3 
for the sub-domain 0 < x < 3. Note that, to limit the 
magnitude of the velocity, the external force is non-zero 
only when the maximum velocity in the domain is smaller 
than 0.3c. For the numerical simulation a = 100, F = 
4/3, 5x/5t = 2.5-\/2, r = 1.0 and a = 0.1 are considered. 

The results of the 3D simulation are shown in FigfTT] 
One can see the initial condition of the magnetic field 
lines in Fig llll a). The colors on the outer domain bound¬ 
aries indicate the value of the vorticity, which is zero ev¬ 
erywhere in the beginning. After applying the external 


We have developed a relativistic MHD lattice Boltz¬ 
mann model, capable of dealing with problems in the 
resistive and ideal regimes. The model is based on the 
relativistic LB model proposed in Ref.jl^, to solve the 
hyd rody namics equations, and the model proposed in 
Ref.[^ to solve the Maxwell equations, where several 
modifications and extensions are implemented to couple 
the models and to use them in the relativistic MHD con¬ 
text. Thus, a D3Q19 lattice configuration is used for the 
hydrodynamic part and a D3Q13 lattice for the electro¬ 
magnetic part. 

The numerical method is validated for test simulations 
in two different regimes, namely propagation of an Alfven 
wave in the ideal MHD limit (high conductivity), and 
evolution of a current sheet in a resistive regime (low con¬ 
ductivity) . The results are compared with the analytical 
ones and very good agreement is observed. Additionally, 
the magnetic reconnection driven by the relativistic KH 
instability is studied in detail and the effect of different 
parameters on the reconnection rate is investigated. It 
is concluded that, while the density ratio has negligible 
effects on the reconnection rate, an increase in the value 
of the shear velocity will decrease the reconnection rate. 
We have also found that, the reconnection rate is pro¬ 
portional to which agrees with the scaling law of 

Sweet-Parker model. Finally, we have presented the re¬ 
sults of 3D simulation of the magnetic reconnection in 
a stellar flare, which is driven by shear velocity in the 
photosphere. We have shown that due to the shear ve¬ 
locity the reconnection happens between flux rope and 
background magnetic field lines. 

It is worth mentioning that, despite the fact that the 
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model in Ref. [2^ is numerically robust at high veloc¬ 
ities, the current model becomes numerically unstable 
when the fluid velocity is high, i.e., in the relativistic su¬ 
personic regime. The reason is not fully clear to us at 
the moment and the issue can be an interesting subject 
for future extension of the work. Also, not that the rel¬ 
ativistic lattice Boltzmann models show nearly an order 
of magnitude faster performance than the corresponding 
hydrodynamic codes [2l|, [2^. Furthermore, the lattice 
Boltzmann model to solve the Maxwell equations also 
turns out to be very efficient, almost an order of magni¬ 
tude faster than Yee’s original FDTD method [IJ. Thus, 
we expect our model to have a competitive performance 


comparing with current models for resistive relativistic 
MHD. 


ACKNOWLEDGMENTS 

We acknowledge financial support from the Euro¬ 
pean Research Council (ERG) Advanced Grant 319968- 
FlowCCS and financial support of the Eidgenssische 
Technische Hochschule Zurich (ETHZ) under Grant No. 
0611-1. 


[1] I. F. Mirabel and L. F. Rodrguez, 
Annu. Rev. Astron. Astrophys. 37, 409 (1999) 

[2] R. Antonucci, Annu. Rev. Astron. Astrophys. 31, 473 
(1993). 

[3] T. Piran, Rev. Mod. Phys. 76, 1143 (2005) 

[4] N. F. Camus, S. S. Komissarov, N. Bucciantini, and P. A. 
Hughes, Mon. Not. R. Astron. Soc. 400, 1241 (2009) 

[5] Y. Masada, S. Nagataki, K. Shibata, and T. Terasawa, 
Publ. Astron. Soc. Japan 62, 1093 (2010). 

[6] S. Komissarov, Mon. Not. R. Astron. Soc. 308, 1069 
(1999). 

[7] M. D. Duez, Y. T. Liu, S. L. Shapiro, and B. C. Stephens, 
Phys. Rev. D 72, 024028 (2005). 

[8] B. D. Farris, T. K. Li, Y. T. Liu, and S. L. Shapiro, 
Phys. Rev. D 78, 024023 (2008). 

[9] J. A. Faber and F. A. Rasio, Living Rev. Relat. 15 (2012). 

[10] A. MacFadyen and S. Woosley, Astrophys. J. 524, 262 
(1999). 

[11] C. Jaroschek, H. Lesch, and R. Treumann, Astrophys. 
J. Lett. 605, L9 (2004). 

[12] J. G. Kirk, O. Skjraasen, and Y. A. Gallant, Astron. 
Astrophys. 388, L29 (2002). 

[13] G. Drenkhahn, Astron. Astrophys. 387, 714 (2002). 

[14] S. Zenitani, M. Hesse, and A. Klimas, Astrophy. J. 696, 
1385 (2009). 

[15] S. S. Komissarov, Mon. Not. R. Astron. Soc. 382, 995 
(2007). 

[16] N. Watanabe and T. Yokoyama, Astrophys. J. Lett. 647, 
L123 (2006). 

[17] C. Palenzuela, L. Lehner, O. Reula, and L. Rezzolla, 
Mon. Not. R. Astron. Soc. 394, 1727 (2009). 

[18] M. Dumbser and O. Zanotti, J. Comput. Phys. 228, 6991 
(2009). 

[19] M. Takamoto and T. Inoue, Astrophy. J. 735, 113 (2011). 

[20] Y. Mizuno, Astrophys. J. Suppl. Ser. 205, 7 (2013). 

[21] M. Mendoza, B. Boghosian, H. Herrmann, and S. Succi, 
Phys. Rev. Lett. 105, 014502 (2010). 

[22] M. Mendoza, I. Karlin, S. Succi, and H. Herrmann, Phys. 
Rev. D 87, 065027 (2013). 

[23] F. Mohseni, M. Mendoza, S. Succi, and H. Herrmann, 
Phys. Rev. D 87, 083003 (2013). 

[24] M. Mendoza and J. Munoz, Phys. Rev. E 82, 056708 

( 2010 ). 

[25] S. Succi, M. Vergassola, and R. Benzi, Phys. Rev. A 43, 
4521 (1991). 


[26] P. J. Dellar, J. Comput. Phys. 179, 95 (2002). 

[27] G. Breyiannis and D. Valougeorgis, Phys. Rev. E 69, 
065702 (2004). 

[28] S. Chen, H. Chen, D. Martnez, and W. Matthaeus, Phys. 
Rev. Lett. 67, 3776 (1991). 

[29] D. O. Martinez, S. Chen, and W. H. Matthaeus, Phys. 
Plasmas 1, 1850 (1994). 

[30] W. Kelvin, Mathematical and Physical Papers Vol. 
4 (Cambridge University Press, Cambridge, England, 
1910) pp. 76-85. 

[31] H. Helmholtz, Philos. Mag., Ser.4, 36, 337 (1868). 

[32] M. Faganello, F. Califano, F. Pegoraro, T. Andreussi, 
and S. Benkadda, Plasma Phys. Contr. F. 54, 124037 
( 2012 ). 

[33] A. Lobanov and J. Zensus, Science 294, 128 (2001). 

[34] M. Vietri, A. Ferrara, and F. Miniati, Astrophy. J. 483, 
262 (1997). 

[35] C.-Y. Wang and R. A. Chevalier, Astrophys. J. 549, 1119 
( 2001 ). 

[36] R. Keppens, G. Toth, R. Westermann, and J. Goed- 
bloed, J. Plasma Phys. 61, 1 (1999). 

[37] T. Nakamura and M. Fujimoto, Adv. Space Res. 37, 522 
(2006). 

[38] H. R. Takahashi, T. Kudoh, Y. Masada, and J. Mat- 
sumoto, Astrophys. J. Lett. 739, L53 (2011). 

[39] C. Cercignani and G. Kremer, The Relativistic Boltz¬ 
mann Equation: Theory and Apllications (Birkhauser, 
Boston; Basel; Berlin, 2002). 

[40] S. Qamar and G. Warnecke, J. Comput. Phys. 205, 182 
(2005). 

[41] D. Ryu, I. Chattopadhyay, and E. Choi, Astrophys. J. 
Suppl. Ser. 166, 410 (2006). 

[42] J. Anderson and H. Witting, Physica 74, 466 (1974). 

[43] D. Hupp, M. Mendoza, I. Bouras, S. Succi, and H. Her¬ 
rmann, Phys. Rev. D 84, 125015 (2011). 

[44] P. A. Sweet, in Electromagnetic phenomena in cosmical 
physics, Vol. 6 (1958) p. 123. 

[45] E. N. Parker, J. Geophys. Res. 62, 509 (1957). 

[46] H. E. Petschek, Magnetic field annihilation, Tech. Rep. 
(DTIC Document, 1963). 

[47] Y. Lyubarsky, Mon. Not. R. Astron. Soc. 358, 113 
(2005). 

[48] S. Zenitani, M. Hesse, and A. Klimas, Astrophys. J. Lett. 
716, L214 (2010). 

[49] E. A. Johnson, Gaussian-Moment Relaxation Closures 



16 


for Verifiable Numerical Simulation of Fast Magnetic Re¬ 
connection in Plasma, Ph.D. thesis, University of Wis¬ 
consin (2013). 

[50] Q. Chen, A. Otto, and L. Lee, J. Geophys. Res. A: Space 


Phys. 102, 151 (1997). 

[51] P. Chen and K. Shibata, Astrophy. J. 545, 524 (2000) 


